program epp;
uses common,ppot;
Const allatoms = 5000;
Type
   v3	    = array[1..3] of myreal;
   v3x3	    = array[1..3] of v3;
   i3	    = array[1..3] of integer;
   r_xyz    = array[1..allatoms] of v3;
   r_e	    = array[1..allatoms] of myreal;
   r_tp	    = array[1..allatoms] of integer;
Var
   i,j,k,l,o,dbg	 : integer;
   txt			 : text;
   r			 : r_xyz;
   ene,ela		 : r_e;
   bl,blt		 : v3x3;
   itp,ito		 : r_tp;
   nat,nats,nsteps,omax	 : integer;
   ich			 : array[1..110] of integer;
   xgamma,dr2,drmax,elas : myreal;
   elaso,epoto		 : myreal;
   dostep		 : boolean;
   uo			 : array[1..100,1..110] of myreal;
   no			 : array[1..100,1..110] of integer;

procedure mxle3(A:v3x3;B:v3;Var v:v3);
Var o,i,j:integer;
    det,dx1,dx2,dx3:real;
begin
   v[1]:=0; v[2]:=0; v[3]:=0;
   DET:=A[1,1]*A[2,2]*A[3,3]+A[1,2]*A[2,3]*A[3,1]
   +A[2,1]*A[3,2]*A[1,3]-A[1,3]*A[2,2]*A[3,1]
   -A[1,1]*A[2,3]*A[3,2]-A[3,3]*A[1,2]*A[2,1];
   if abs(det)<0.000001 then begin
      writeln('MXLE3 error: det=0'); halt;
   end;
   DX1:=B[1]  *A[2,2]*A[3,3]+A[1,2]*A[2,3]*B[3]
   +B[2]  *A[3,2]*A[1,3]-A[1,3]*A[2,2]*B[3]
   -B[1]  *A[2,3]*A[3,2]-A[3,3]*A[1,2]*B[2];
   DX2:=A[1,1]*B[2]  *A[3,3]+B[1]  *A[2,3]*A[3,1]
   +A[2,1]*B[3]  *A[1,3]-A[1,3]*B[2]  *A[3,1]
   -A[1,1]*A[2,3]*B[3]  -A[3,3]*B[1]  *A[2,1];
   DX3:=A[1,1]*A[2,2]*B[3]  +A[1,2]*B[2]  *A[3,1]
   +A[2,1]*A[3,2]*B[1]  -B[1]  *A[2,2]*A[3,1]
   -A[1,1]*B[2]  *A[3,2]-B[3]  *A[1,2]*A[2,1];
   v[1]:=dx1/det; v[2]:=dx2/det; v[3]:=dx3/det;
end;
function bound(z: myreal):myreal;
begin
   bound:=z;
   if z>0.4999999 then bound:=z-1.0;
   if z<-0.5000001 then bound:=z+1.0;
end; { bound }
procedure relax(move :boolean; nat : integer; itp:r_tp; xgamma,drmax:myreal;Var r:r_xyz; Var e:r_e; Var epot,elas,dr2 : myreal);
Var
   i,j,k,o,m,n	      : integer;
   bbij,ar2,cf	      : myreal;
   drt		      : v3;
   bb		      : array[1..allatoms] of v3x3;
   io,jo	      : integer;
   i1,i2,p1,p2,p3,ina : integer;
   p		      : i3;
   ar,u,ud,udd,drm2   : myreal;
   v,w		      : v3;
   bo		      : boolean;
   b,dra	      : r_xyz;
begin
   drm2:=drmax*drmax;
   epot:=0.0; ina:=0; dr2:=0.0; elas:=0;
   for i:=1 to nat do begin
      l:=0; ene[i]:=0.0; if move then ela[i]:=0.0;
      for io:=1 to 3 do begin
	 b[i,io]:=0.0;
	 for jo:=1 to 3 do bb[i,io,jo]:=0.0;
      end;
   end;
   for i:=1 to nat do begin
      for j:=i to nat do begin
	 for p1:=-iext[1] to iext[1] do
	    for p2:=-iext[2] to iext[2] do
	       for p3:=-iext[3] to iext[3] do begin
		  bo:=((p1=0) and (p2=0) and (p3=0) and (i=j));
		  if not(bo) then begin
		     for o:=1 to 3 do v[o]:=bound(r[j,o]-r[i,o]);
		     p[1]:=p1; p[2]:=p2; p[3]:=p3;
		     for m:=1 to 3 do w[m]:=0.0;
		     for m:=1 to 3 do
			for n:=1 to 3 do w[n]:=w[n]+(v[m]+p[m])*bl[m,n];
		     ar2:=w[1]*w[1]+w[2]*w[2]+w[3]*w[3];
		     ar:=sqrt(ar2);
		     if (ar<rcut) and (ar>0.01) then begin
			l:=l+1;
			upp(itp[i],itp[j],ar,u,ud,udd);
			epot:=epot+u;
			if (i=j) then epot:=epot-u/2.0;
			ene[i]:=ene[i]+u/2.0;
			ene[j]:=ene[j]+u/2.0;
			for io:=1 to 3 do begin
			   b[i,io]:=b[i,io]+ud*w[io];
			   for jo:=1 to 3 do begin
			      bbij:=(ud-udd)*w[io]*w[jo]/ar2;
			      IF (IO=JO) then BBIJ:=BBIJ-ud;
			      BB[i,IO,JO]:=BB[i,IO,JO]+BBIJ;
			   end;
			end;
			for o:=1 to 3 do w[o]:=-w[o];
			for io:=1 to 3 do begin
			   b[j,io]:=b[j,io]+ud*w[io];
			   for jo:=1 to 3 do begin
			      bbij:=(ud-udd)*w[io]*w[jo]/ar2;
			      IF (IO=JO) then BBIJ:=BBIJ-ud;
			      BB[j,IO,JO]:=BB[j,IO,JO]+BBIJ;
			   end;
			end;
		     end;
		  end;
	       end;
      end;
      if l>ina then ina:=l;
   end;
   for i:=1 to nat do begin
      mxle3(bb[i],b[i],drt);
      for io:=1 to 3 do
	 for jo:=1 to 3 do
	    ela[i]:=ela[i]+drt[io]*bb[i,io,jo]*drt[jo];
      ela[i]:=-ela[i]/2;
      elas:=elas+ela[i];
      if move then begin
	 ar2:=0.0; for o:=1 to 3 do ar2:=ar2+drt[o]*drt[o];
	 cf:=1.0;
	 if ar2>drm2 then begin
	    cf:=sqrt(drm2/ar2);
	    for o:=1 to 3 do drt[o]:=drt[o]*cf;
	 end;
	 dr2:=dr2+ar2*cf*cf;
	 for o:=1 to 3 do drt[o]:=drt[o]*xgamma;
	 mxle3(blt,drt,dra[i]);
	 for o:=1 to 3 do r[i,o]:=bound(r[i,o]-dra[i,o]);
      end;
   end;
end; { relax }
begin
   readln(nsteps,xgamma,drmax);
   dbg:=0;
   pp_inp;
   for i:=1 to 3 do begin
      for o:=1 to 3 do read(bl[i,o]); readln;
      for o:=1 to 3 do blt[o,i]:=bl[i,o];
   end;
   omax:=0;
   readln(nats);
   for i:=1 to nats do begin
      nat:=nat+1;
      for o:=1 to 3 do read(r[nat,o]); readln(itp[nat],ito[nat]);
      if ito[nat]>omax then omax:=ito[nat];
      if itp[nat]=0 then nat:=nat-1;
   end;
   for i:=1 to 110 do ich[i]:=0;
   for i:=1 to nat do ich[itp[i]]:=ich[itp[i]]+1;
   for i:=1 to 100 do for j:=1 to 110 do begin  {orbit energies}
      uo[i,j]:=0.0; no[i,j]:=0;
   end;
   writeln('step E[eV/at] E_ela[eV/at] dE dE_ela dR[A] ');
   dostep:=true;
   epoto:=0.0; elaso:=0.0;
   for i:=1 to nsteps do begin
      relax(dostep,nat,itp,xgamma,drmax,r,ene,epot,elas,dr2);
      if(i=1) then begin
	 writeln((i-1):1,' ',epot/nat:1:8,' ',elas/nat:1:10,'  --  -- ',sqrt(dr2/nat):1:10);
	 elaso:=elas;
	 epoto:=epot;
      end else
	 writeln((i-1):1,' ',epot/nat:1:8,' ',elas/nat:1:10,' ',(epot-epoto)/nat:1:8,' ',(elas-elaso)/nat:1:8,' ',sqrt(dr2/nat):1:10);
   end;
   dostep:=false;
   relax(dostep,nat,itp,xgamma,drmax,r,ene,epot,elas,dr2);
   writeln((nsteps):5,' ',epot/nat:15:8,' ',sqrt(dr2/nat):10:5);
   {writeln('E_rel=',epot/nat:1:8,'[eV/atom]');}
   rewrite(txt,'xyz.epp');
   for i:=1 to 3 do begin
      for j:=1 to 3 do write(txt,bl[i,j]:12:5); writeln(txt);
   end;
   writeln(txt,nat:1,' # atoms (no vacancies)');
   for i:=1 to nat do begin
      for o:=1 to 3 do write(txt,r[i,o]:1:8,' ');
      writeln(txt,'   ',itp[i]:2,' ',ito[i]:2,' ',ene[i]:1:8,' ',ela[i]:1:8);
   end;
   close(txt);
   {orbit energies}
   for i:=1 to nat do begin
      uo[ito[i],itp[i]]:=uo[ito[i],itp[i]]+ene[i];
      no[ito[i],itp[i]]:=no[ito[i],itp[i]]+1;
   end;
   {for i:=1 to omax do for j:=1 to 110 do if no[i,j]>0 then writeln(i:2,' ',j:2,' ',uo[i,j]/no[i,j]:1:4);}
end.
